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The temperature-dependent Matsubara Green’s function that is used to describe temperature- 
dependent behavior is expressed on a numerical grid. While such a grid usually has a couple of 
hundred points for low-energy model systems, for realistic systems in large basis sets the size of an 
accurate grid can be tens of thousands of points, constituting a severe computational and memory 
bottleneck. In this paper, we determine efficient imaginary time grids for the temperature-dependent 
Matsubara Green’s function formalism that can be used for calculations on realistic systems. We 
show that due to the use of orthogonal polynomial transform, we can restrict the imaginary time 
grid to few hundred points and reach micro-Hartree accuracy in the electronic energy evaluation. 

Moreover, we show that only a limited number of orthogonal polynomial expansion coefficients are 
necessary to preserve accuracy when working with a dual representation of Green’s function or 
self-energy and transforming between the imaginary time and Matsubara frequency domain. 


I. INTRODUCTION 

The use of numerical grids in calculations for realistic 
systems has a long history in quantum chemistry simula¬ 
tions. For example, in density functional theory (DFT) 
a numerical integration is necessary for the evaluation 
of the exchange-correlation contribution to the density 
functionab^®. Similarly, in Laplace transformed MP2 
(LT-MP2) a quadrature is used to represent an integral 
that leads to the removal of the energy denominators®^^^. 
Recently, an implementation of random phase approx¬ 
imation (RPA)^^ appeared that uses an efficient imag¬ 
inary time grid to yield temperature-independent RPA 
energy. The above mentioned methods are just a few 
examples of using efficient quadrature, a more exten¬ 
sive literature on the subject can be found in refs 13- 
32. Thus, it is fair to say that extensive knowledge exist 
on representing temperature-independent quantities on 
a grid when ground state methods are used. However, 
very little is known about how to efficiently represent 
temperature-dependent data on finite-temperature imag¬ 
inary axis Matsubara grids. 

Several factors distinguish the finite temperature 
Green’s function from the zero-temperature Green’s func¬ 
tion formalism. Firstly, let us note that the temperature- 
dependent Green’s function is a discrete object for which 
the grid points itUn are spaced according to the Mat¬ 
subara grid Wn = {2n + 1 )tt/(3, where n G Z and 
/3 = l/k^T is the inverse temperature. In compari¬ 
son, the zero-temperature Green’s function represented 
on imaginary axis is a continuous function. Similarly, the 
temperature-dependent imaginary time Green’s function 
is an antiperiodic function between 0 and /3, while the 
zero-temperature Green’s function of imaginary time is 
a non-periodic function decaying rapidly and smoothly 
to zero. Consequently, traditional quadratures devel¬ 
oped for zero-temperature RPA Green’s functions used 
in ref 12 or LT-MP2 to represent denominators® are not 
suitable for temperature-dependent Green’s function cal¬ 
culations. 


Currently, temperature-dependent Green’s function 
calculations are mostly done for low-energy model sys¬ 
tems. While such calculations for large realistic systems 
are still in their infancy, one can easily imagine that they 
could be very important in materials science for mate¬ 
rials with small band gaps where the change of proper¬ 
ties with temperature is significant and multiple states 
can be easily populated, or for a system which exhibit 
temperature-dependent phase transition caused by the 
electronic degrees of freedom. 

Unlike the low-energy models, the orbital energies in 
realistic systems span a huge energy window frequently 
varying even between —3000 eV to 300 eV, see Figure 1. 
Thus, when quantitative accuracy in the calculations of 
realistic systems is desired, new challenges arise that are 
not present in model system calculations. 

The temperature-dependent Green’s functions have to 
be represented on a numerical grid with a spacing de¬ 
fined by the temperature but covering the energy win¬ 
dow spanned by the realistic system. This usually results 
in a grid containing thousands of frequencies. Conse¬ 
quently, the Fourier transform from the imaginary time 
to the imaginary frequency axis, according to Nyquist 
theorem, requires twice as many imaginary time points as 
frequency points to yield accurate results. Having several 
thousands data points makes the calculations for realistic 
systems extremely challenging even if each grid point can 
be calculated in parallel. 

Motivated by the aforementioned challenges, we first 
determine how the imaginary time grid can be truncated 
to a reasonable size and how the resulting errors can 
be controlled. To achieve it, we replace the numerical 
Fourier transform (ir -G iuj) with an orthogonal poly¬ 
nomial transform (ir -G L) and {L -G ioj), where as L 
we denote expansion coefficients of a Green’s function 
or self-energy in an orthogonal polynomial basis. Sub¬ 
sequently, we examine if (i) the expansion coefficients L 
can be produced using a smaller number of grid points 
than currently used in the numerical Fourier transform, 
and if (ii) the number of expansion coefficients L is small 
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FIG. 1. Orbital energies for atoms and molecules in large 
basis sets containing diffuse orbitals. 


enough that they can be easily stored in realistic calcu¬ 
lations. 

This paper is organized as follows. In section II, we 
review the necessary theoretical background leading to 
compact imaginary time grids and orthogonal polynomial 
transforms (tr —>■ L —?► ico) using Legendre expansion co¬ 
efficients. In section III, we discuss how an orthogonal 
polynomial transform can be used in the second order 
iterative Green’s function method (GF2) to reduce the 
size of the imaginary time grid. In section IV, we present 
numerical results showing that the required number of 
imaginary time points is much smaller than in the origi¬ 
nal uniform and power-law grid and that the number of 
required expansion coefficients can be kept small even if 
the micro-Hartree (/rHa) accuracy is desired. Finally, we 
present conclusions in section V. 


In many Green’s function methods two Green’s func¬ 
tion (G) or self-energy (S) representations are used: an 
imaginary time representation G{iT) or S(fr) and imagi¬ 
nary frequency representation G{i(jj) or Yi(iLo). In an effi¬ 
cient implementation, one frequently changes from ir to 
iui and back, depending if it is more computationally ad¬ 
vantageous to work with a it or iuj representation. Note 
that while the above statement is general and the it to 
iio transform may be present in temperature-independent 
calculations such as RPA or in this paper we 

focus exclusively on the imaginary time and frequency 
used for temperature-dependent Green’s functions^®^^®. 

Thus, frequently a Green’s function method proceeds 
according to scheme illustrated in Figure 2 and the com¬ 
putational bottleneck lies in the evaluation of T,[iT). 

Dyson equation 

G{iuj) =[[GoM]-^ > G(tr) 

Self-energy 

evaluation 

IFFT 

S(za;) - - 

imaginary Matsubara imaginary 

frequency time 

FIG. 2. An example of the representation change between 
imaginary time and frequency in Green’s functions methods. 
We denote the non-interacting and interacting/correlated 
Green’s functions as Goiiu)) and G{iui), respectively. 


The number of imaginary time points used to represent 
E(fT) enters as a prefactor. However, a sizable prefactor 
makes such calculations for large system impossible. 

There are two reasons for the large size of imaginary 
time grids for temperature-dependent Green’s functions, 
firstly enough points are required to preserve the tem¬ 
perature dependence, secondly a significant number of 
points is required to make a numerical Fourier transform 
accurate. 


A. Orthogonal polynomial representation of 
self-energy 

We consider a single-particle temperature-dependent 
imaginary time self-energy S(it) on the interval [0,/3]. 
The self-energy can be expanded in Legendre-Fourier se¬ 
ries using orthogonal polynomials P/. Note that different 
orthogonal polynomials, e.g. Ghebyshev or Legendre, can 
be used as a basis. However, all results presented here 
should be valid regardless of the representation®®’^®. In 
this work, we use Legendre polynomials defined on the 
interval [—1,1]. The Legendre polynomial expansion of 
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the self-energy on this interval is given by 


B. Sparse imaginary time grid 


/07 _|_ 1 

^- Pi{x{iT))I]\^, ( 1 ) 

Z>0 ^ 

where Pi{x{iT)) is a Legendre polynomial of rank I and 
S-j is the corresponding expansion coefficient, 

= v/2l + T / dTP,(a:(T))Sy(zr), (2) 

Jo 

and x{t) = 2t//3 — 1 maps interval [0, /3] onto [—1,1]- For 
detailed description of the properties of Legendre poly¬ 
nomials see ref 41. 

In the past, several advantages of using an orthogo¬ 
nal polynomial representation of Green’s function (and 
self-energy) were explored. Boehnke et. al.^® used the 
Legendre polynomial representation of Green’s functions 
as a noise filter in continuous-time hybridization expan¬ 
sion quantum Monte Carlo (CT-HYB) calculations. Ad¬ 
ditionally, the Legendre representation of Green’s func¬ 
tion can be more compact than Matsubara Green’s func- 

tion39.42.43^ 

In this paper, we describe another advantage of us¬ 
ing Legendre representation. It allows us to perform an 
orthogonal polynomial two step transform (S(zt) —>■ L) 
and (L Y]{iuj)) and to reduce the size of the imaginary 
time grid that would be otherwise necessary to maintain 
a high accuracy. This (S(ir) -f-)- L) transformation is 
unitary and can be written as a matrix multiplication 


In this section, we examine the number of imaginary 
time grid points necessary to perform the Fourier trans¬ 
form (S(ir) —>■ L) accurately. This assessment is ab¬ 
solutely vital to the success of many approaches where 
evaluation of the self-energy on the time grid is the com¬ 
putation bottleneck. 

To estimate the number of grid points necessary we 
performed calculations for atoms and molecules using dif¬ 
ferent numbers of imaginary time points. While both 
atoms and molecules that we use here as test examples 
do not display different physics for a large range of tem¬ 
peratures due to the size of the gap present in these sys¬ 
tems, they are very challenging examples since we strive 
to calculate electronic energy for very low temperatures 
near absolute zero (large value of inverse temperature 
/3 = 100). Consequently, our grid spacing has to be very 
small and the grid has to span significant energy win¬ 
dow shown in Figure 1 thus requiring very many points 
if the uniform or power-law grids are used. For numeri¬ 
cal Fourier transforms we used two different grids, a uni¬ 
form grid and power-law grid as described below. For the 
(S(jt) —>■ L) transform, we used a modification of the 
power-law grid to compute a fixed number of Legendre 
coefficients. We have chosen 200 expansion coefficients 
since such a number of coefficients allows us to calculate 
the energy to pHa accuracy. 

1. Uniform imaginary time grid 


^ /o/ I 1 rP 

dTP--^Pdx{r)) 


l>0 

oo 


/5 

= YTnlAj, ( 3 ) 

/>o 

where T^i is the unitary matrix with elements defined as 
Tni = , (4) 


where ji{z) are the spherical Bessel functions of the sec¬ 
ond kind. 

Even though the Legendre series is infinite in princi¬ 
ple, in all practical calculations only a finite number of 
expansion coefficients is used. For all atomic and molecu¬ 
lar systems studied here the expansion coefficients decay 
very fast. In the worst case, only a few hundreds of them 
are necessary (see section IV), therefore only 0{Niv?) 
double precision numbers have to be stored, where n is 
the number of orbitals and Ni is the number of terms in 
the Legendre expansion. In contrast to calculations em¬ 
ploying orthogonal polynomial transform, a typical nu¬ 
merical Fourier transform may require tens to hundreds 
of thousands of imaginary time grid points N^, making 
the cost of evaluation of the self-energy very significant. 


The grid points are uniformly spaced within the inter¬ 
val [0,/3] where the antiperiodic G(it) or S(iT) is repre¬ 
sented. 


2. Power-law time grid 

Since imaginary time self-energy is sharply peaked 
around endpoints (0 and ±/3) it is convenient to use 
non-uniformly spaced grids to represent it. A power-law 
grid44 is constructed to be dense around endpoints and 
sparse between them where imaginary time self-energy 
is close to zero. The power-law grid is defined by two 
parameters: the power coefficient p and the uniform co¬ 
efficient u. The first step in creating such a grid is placing 
points with the coordinates r, = fd/2^,j G {0, ...,p— 1} 
starting from each endpoint and also placing a midpoint 
at f3/2. Gonsequently, such a grid has 2p -|- 1 points. 
Then, each interval between power points is divided into 
2u uniformly spaced subintervals. 


3. fS(ir) L) transform with a power-law grid 

We evaluated Legendre expansion coefficients employ¬ 
ing only a fraction of the original power-law grid. We 
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kept number of power points fixed p = 12 and choose the 
number of uniform points as u = 2^ where n € {0,5}. 
The resulting grid has 2u(p+l) + l imaginary time points 
that correspond to 53, 105, 209, 417 and 833 points for 
n G {0,..., 5} respectively. 

Using the uniform, power-law grid and orthogonal 
polynomial transform, we evaluated the Matsubara self¬ 
energy S(iu>). In Figure 3, we show convergence for dif¬ 
ferent grids. We set as our reference the self-energy ob¬ 
tained with 200 Legendre polynomials. 



Log 


FIG. 3. Convergence of the self-energy as a function 

of number of imaginary time points. 


As seen in Figure 3, the (S(ir) — L) and (L -G 
S(ja;)) transforms converge much faster than the regu¬ 
lar (S(iT) —>■ T,(iu})) Fourier transform performed using 
a uniform or power-law grid. 


III. SELF-CONSISTENT SECOND-ORDER 
GREEN’S FUNCTION THEORY USING 
ORTHOGONAL POLYNOMIAL TRANSFORM 


— o+<^+ 




FIG. 4. GF2 is the second order approximation to the self¬ 
energy. It includes two first order self-energy diagrams (from 
left to right): Hartree (direct), Fock (exchange) and two sec¬ 
ond order self-energy diagrams: direct and exchange. 


In this section, as an application to realistic calcu¬ 
lations, we briefly describe the framework of the self- 


consistent second-order Green’s function theory (GF2) 
using an orthogonal polynomial transform. 

GF2 employs the second-order approximation to the 
self-energy resulting in inclusion of all Feynman self¬ 
energy diagrams up to the second order, as shown in 
Figure 4. The first two diagrams are already included 
at mean-field level. Within spin-restricted GF2, the last 
two diagrams are translated into the following expression 
for the imaginary time self-energy 

U*j(Ir„) = - Gli(iTn)Gl^^{iTn)Glq{-iTn) 

klmnpq 

^ '^imqk {^'k^lpnj Vnplj) , (5) 

where G*^(jt„) is the zeroth order imaginary time Green’s 
function and Vijki are two-electron integrals defined as 

Vijki= J J ( 6 ) 

In GF2, as illustrated in Figure 2 for a general case, the 
calculation of the second-order self-energy according to 
eq 5 is done in the imaginary time domain while the 
Dyson equation from Figure 2 is much easier to solve in 
the frequency domain. 

In our previous implementation, a typical molecular 
all electron GF2 calculation that involves both core and 
virtual orbitals may require many thousands of Matsub¬ 
ara frequencies making the total amount of storage 
necessary equal to where n is the number of 

orbitals. The imaginary time Green’s function is repre¬ 
sented by 0{NT'n?) double precision numbers, where 
is the number of points of the imaginary time grid. Build¬ 
ing the self-energy according to eq 5 scales as O(Wn^), 
and despite that the self-energy calculation at any given 
imaginary time point is independent and can be made 
parallel, a large prefactor W is slowing down calcula¬ 
tions significantly even when using a power-law grid. 

Since employing the orthogonal polynomial transform 
restricts the imaginary time grid even for the most dif¬ 
ficult cases to fewer than 400 points, we implemented 
it as part of our algorithm. Here, we give a complete 
step-by-step modified algorithm description. 

1. Start with Hartree-Fock (HF) reference solution 
(although starting from DFT reference is equally 
possible and advantageous for cases that are diffi¬ 
cult to converge using HF) and build initial Mat¬ 
subara Green’s function in non-orthogonal AO ba¬ 
sis according to: 

Go(Ia;„) = [(M + *Wn)S-F]-\ (7) 

where S is the overlap matrix, F is the Fock matrix 
and /r is the chemical potential. 

2. Perform discrete Fourier transform of Go(ja;„) into 
its imaginary time counterpart Go(*'r). Alterna¬ 
tively, at this point, it is possible to avoid discrete 
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Fourier transform if one starts directly from the 
imaginary time HF Green’s function constructed 
from HF orbital energies 

G° (iT„) = 0(*r„) {n{E,) - 1) 

( 8 ) 

where 0{x) is Heaviside step function, Ei = ei — 
and n{Ei) = l/(e^'®‘ +1) is the Fermi distribution. 
Since the Green’s function from eq 8 is constructed 
using MO orbital energies, it should be transformed 
to AO basis before proceeding to the next step. 

3. Galculate the self-energy on the imaginary time 
grid according to eq 5. It is at this point where we 
first take advantage of Legendre polynomial rep¬ 
resentation of self-energy, since the Legendre rep¬ 
resentation allows us to use small imaginary time 
grids with only a fraction of points of the grid we 
used in our original implementation^®. 

4. Obtain the Legendre expansion coefficients by per¬ 
forming an integration of the self-energy S(it) ac¬ 
cording to eq 2. 

5. Build the imaginary frequency self-energy S(*a;) by 
performing a transform of Legendre coefficients ac¬ 
cording to eq 4. 

6. Solve the Dyson equation to obtain an updated 
Green’s function. 

7. Find the chemical potential fj, to ensure that a 
proper number of electrons is present in the sys¬ 
tem. 

8. Galculate the density matrix and use it to update 
Fock matrix. 

9. Go to point 6 and iterate until the density matrix 
and chemical potential /i stop to change. 

10. Galculate the one-body energy as 

Eib = —- + Fjuy), (9) 

where = —2Gu^{iT = (3) is the correlated den¬ 
sity matrix and the Fock matrix is evaluated using 
this correlated density matrix. The two-body en¬ 
ergy can be evaluated using as 

^ l^freg 

E 2 b = - o X! Tr[G(*a;„)E(iw„)], (10) 

^ n=0 

for details see ref 45. 

11. Transform G{ioj) to G(ir) and go to step 3 and 
iterate until the total energy converges. 


In this section, we provide results of atomic and molec¬ 
ular calculations with the above introduced GF2 algo¬ 
rithm. To assess the accuracy and efficiency of the al¬ 
gorithm described above, we performed several bench¬ 
mark calculations with large basis sets using diffuse or¬ 
bitals. These basis sets usually required the most ex¬ 
tensive imaginary time grid and are necessary to reach 
quantitative accuracy and converge with basis set size. 

Additionally, we also tested our algorithm on a few sys¬ 
tems with transition metal atoms with ecp-sdd-DZ^®^^® 
basis set containing pseudopotentials for inner shell elec¬ 
trons. We investigated systems with pseudopotentials 
since these are frequently used in solid state calculations 
and it is our interest to assess how compact the grids can 
become for such systems. 

Our investigations can be divided into two groups 
evoking our original questions about (i) the size of the 
grid necessary to calculate the Legendre coefficients ac¬ 
curately and (ii) the compactness of the Legendre ex¬ 
pansion. 



FIG. 5. Convergence of electron correlation energy with the 
size of imaginary time grid Nr calculated for a few atoms and 
molecules. 


A. Convergence of the electronic energy with 
respect to the grid size 

We calculated the electronic correlation energy accord¬ 
ing to eq 10 for the grids defined in Section HB3 and 
compared it to the one obtained using our previous GF2 
implementation with sufhciently large imaginary time 
grids. Our previous implementation required at least an 
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FIG. 6. Even coefficients of Legendre expansion of self-energy calculated using imaginary time grid consisting of 417 points for 
atoms and molecules. Left panel shows zoomed in region where E* ~ 0 where Gibbs oscillations start to develop. 


order of magnitude larger grids than the current one us¬ 
ing the Legendre expansion. 

In Figure 5, we plot the error in electronic correlation 
energy obtained by using different number of imaginary 
time points that are used to produce Legendre expansion 
of self-energy consisting of 200 expansion coefficients. Let 
us first note that using 200 Legendre coefficients results 
on average in 20 or less /rHa error in the total energy, see 
Table 1. 

If we set the value of correlation energy using 200 Leg¬ 
endre polynomials as a reference then it is evident from 
Figure 5 that with only 104 imaginary time grid points 
an acceptable accuracy (less than 50 /rHa from the exact 
answer) can be achieved. This is already far below the 
commonly accepted chemical accuracy of « 1 kcal/mol. 
Using fewer than 100 grid points is not advisable since for 
53 points for every examined case the accuracy was about 
0.5 mHa away from the exact answer that is unaccept¬ 
able in almost any quantum chemical calculation. For 
imaginary time grids with 209 and more points the accu¬ 
racy reaches a plateau since for 200 Legendre coefficients 
a grid of 209 frequency points is sufficient to produce 
these coefficients with accuracy reaching numerical pre¬ 
cision. If we desired to reach better accuracy than /iHa 
level, then a larger number of the Legendre expansion 
coefficients should be employed in our calculation. 


B. Compactness of the Legendre expansion 

Next, we will study how the accuracy of our calcula¬ 
tions depends on the number of terms in Legendre ex¬ 
pansion of the self-energy. We fix the size of imaginary 
time grid to 417 points because, as it is seen from Fig¬ 


TABLE 1. Error in /iHa of the total energy for different num¬ 
ber of Legendre expansion coefficients for several atoms and 
molecules. 


atom or molecule 

1=200 

1=100 

1=40 

He 

0.09 

0.05 

-44.47 

Be 

0.30 

0.02 

-8.29 

Ne 

-1.05 

-69.57 

-177.23 

Mg 

0.50 

-6.70 

-6.28 

Ar 

-9.93 

-28.53 

697.07 

Ca 

-3.42 

-7.54 

401.17 

H20“ 

-22.00 

-99.99 

9478.18 

(H20)2“ 

-27.26 

372.33 

21048.60 

(H20)3“ 

-53.20 

47.81 

29737.13 

(H20)4“ 

-71.21 

63.47 

39668.95 

NHa* 

-7.12 

-225.17 

6953.89 

CeHs" 

-7.87 

-668.08 

28778.84 

Cd 

-2.27 

-38.27 

201.04 

CU2“ 

-31.09 

-151.26 

6794.16 

Pd 

-2.25 

-48.21 

121.04 

Ag2" 

-4.67 

-36.86 

658.88 


“ Geometry was taken from ref 50. 

* Experimental geometries were taken from NIST 
Gomputational Chemistry Comparison and Benchmark 
Database®'^. 
d(Cu-Cu)=3.63 a. u. 
d(Ag-Ag)=5.46 a. u. 


ure 5, this number of imaginary frequencies Is sufficient 
to produce at least 200 accurate Legendre coefficients in 
the expansion of the self-energy. 

First, we will look at values of expansion coefficients 
Eqq determined by integrating self-energy according to 
eq 2. In Figure 6 we plotted even Legendre coefficients for 
a few atoms and molecules. Odd coefficients show very 
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FIG. 7. Difference between the standard GF2 correlation en¬ 
ergy and correlation energy obtained by using Legendre ex¬ 
pansion of the self-energy for several closed-shell atoms with 
aug-cc-pVDZ basis set (cc-pVDZ basis set was used for Ca). 



FIG. 8. Difference between the standard GF2 correlation en¬ 
ergy and correlation energy obtained by using Legendre ex¬ 
pansion of the self-energy for several transition metal atoms 
and diatomic clusters with ecp-sdd-DZ basis set. 


similar behavior in almost all cases. As seen from the 
left panel of Figure 6, Legendre coefficients decay mono- 
tonically, converging to zero with the decay rate that is 
system-specihc. For all cases studied in this work, we 
observe a fast decay of Legendre expansion coefficients 
used to represent the imaginary time self-energy. Thus, 
the Legendre polynomials form a compact representation 



FIG. 9. Difference between the standard GF2 correlation en¬ 
ergy and correlation energy obtained by using Legendre ex¬ 
pansion of the self-energy for several small molecules with TZ 
basis set. 


not only for Hubbard model^® but also for atoms and 
molecules. Since realistic molecular systems have dif¬ 
fuse orbitals and span large energy spectrum an increase 
in the number of expansion terms in comparison to the 
Hubbard model is to be expected. 

The right panel of Figure 6 shows a zoom in the re¬ 
gion where the expansion coefficients are close to zero, 
Eqq = 0. a closer look at the values of Eqq reveals a 
numerical noise. This noise is known as Gibbs’s oscil¬ 
lations®^ and arises when too few imaginary time grid 
points are used to evaluate higher order Legendre coef¬ 
ficients. In order to prevent numerical noise buildup af¬ 
fecting very high orders of the Legendre expansion, one 
should truncate the expansion once oscillations are de¬ 
tected. Another method is to damp Gibbs oscillations 
by introducing an integral kernel function®^. This option 
has been previously explored in the context of the Hub¬ 
bard model^'^. The particular choice of integral kernel 
function depends on several factors and is not convenient, 
especially if a black-box method is desired. Motivated 
by our aspiration to keep the self-energy as compact as 
possible, we took the first route and studied how the 
truncation of the Legendre expansion of the self-energy 
influences the accuracy of the method. The truncation 
of the expansion is unambiguous and the resulting self¬ 
energy is free from Gibbs oscillations. It is also worth 
mentioning that truncation criteria can be easily imple¬ 
mented and do not require any special care from the user 
and hence can be introduced as a part of any black-box 
computational package. 

We performed calculations for our test set containing 
atoms and simple molecules truncating the Legendre se¬ 
ries after various number of terms = 0 for Z > lent- 










Figures 7-9 show that only couple hundreds of Leg¬ 
endre expansion coefficients are necessary to yield very 
accurate results. All these results were calculated using 
417 imaginary time points. In all cases very fast con¬ 
vergence was achieved and less then 100 Legendre poly¬ 
nomials were needed to converge correlation energy to 1 
ruHa. We observe that for all atoms considered here the 
energy continues to converge to the exact one and about 
200 Legendre polynomials are needed to recover it up to 
the /iHa from the exact result. Similar observations can 
be made for calculations involving pseudopotentials and 
thus our method can be reliably applied to calculations of 
complex systems containing both transition metals and 
light atoms. A somewhat slower convergence was ob¬ 
served in the case of molecular calculations. In this case 
more than 200 but less than 300 Legendre polynomials 
were necessary to get within few /rHa from the exact en¬ 
ergy. 

V. SUMMARY AND CONCLUSIONS 

The frequency—time duality is present in many 
temperature-dependent methods and G'(*t) or S(ir) can 
be transformed to G(iuj) or E(ia;) depending if handling 
the time or frequency object is more computationally ad¬ 
vantageous. Both the frequency and time Green’s func¬ 
tions have to be represented on a numerical grid, thus a 
small number of grid points is crucial for achieving com¬ 
putational efficiency. 

While the construction of imaginary time grids 
is a well-studied problem and appears in Laplace- 
transformed Mpller-Plesset (LT-MP2) perturbation the¬ 
ory, such grids are appropriate for the zero-temperature 
Green’s function but cannot be employed to study the 
temperature-dependent Green’s function. 

In this paper, we have presented a method that makes 
transform between imaginary time and imaginary fre¬ 


quency for temperature-dependent Green’s function con¬ 
verge much faster with respect to the number of neces¬ 
sary imaginary time points than the traditional uniform 
and power-law grid. To achieve this goal we have used 
a combination of a very sparse power-law grid together 
with explicit transform based on a Legendre expansion 
of the self-energy. We have shown that to converge the 
Legendre coefficients necessary to perform the explicit 
transform we need an order of magnitude fewer imagi¬ 
nary time points that when doing the numerical Fourier 
transform. 

Moreover, we have also shown that even for realistic 
systems in basis sets with a significant energy spread only 
a few hundred (200-300) of Legendre coefficients are nec¬ 
essary to reach the accuracy of ^Ha when compared with 
the fully converged result. Overall the orthogonal poly¬ 
nomial representation of self-energy offers fast, accurate, 
and less storage demanding solution for practical realistic 
calculations. 

We have also applied such a representation of the self¬ 
energy to the GF2 method resulting in very accurate en¬ 
ergies for atoms and molecules while using only a limited 
number of imaginary time grid points and only couple 
hundred (200-300) of Legendre expansion coefficients. 

While at present no large scale realistic calculations are 
performed including temperature coming from the elec¬ 
tronic effects, due to the increasing interest in new ma¬ 
terials, we believe that such calculations will become im¬ 
portant in the near future and identifying and overcoming 
major bottlenecks connected to the efficient representa¬ 
tion of Green’s function and self-energy in the imaginary 
time domain is an important step in this new direction. 
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